Practical 6

Ref URLs:

http://openplaques.org/

An Introduction to Analysing Spatial Patterns

1. Point Pattern Analysis

The question we want to answer is: “For any given London Borough, are the Blue Plaques within that borough distributed randomly or do they exhibit some kind of dispersed or clustered pattern?”

libraries

#install.packages("spatstat")
#install.packages("GISTools")
library(spatstat)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
library(tmap)
library(sf)
library(geojsonio)
library(tmaptools)

Set up data

London Boroughs

##First, get the London Borough Boundaries
EW <- geojson_read("http://geoportal.statistics.gov.uk/datasets/8edafbe3276d4b56aec60991cbddda50_2.geojson", what = "sp")
#pull out london using grep and the regex wildcard for'start of the string' (^) to to look for the bit of the district code that relates to London (E09) from the 'lad15cd' column in the data slot of our spatial polygons dataframe
BoroughMap <- EW[grep("^E09",EW@data$lad15cd),]
#plot it using the base plot function
qtm(BoroughMap)

Blue Plaques

##Now get the location of all Blue Plaques in the City
BluePlaques <- geojson_read("https://s3.eu-west-2.amazonaws.com/openplaques/open-plaques-london-2018-04-08.geojson", what = "sp")

Check CRS

summary(BoroughMap)

Object of class SpatialPolygonsDataFrame
Coordinates:
        min        max
x -0.510277  0.3340243
y 51.286759 51.6918756
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
      lad15cd                   lad15nm            lad15nmw     objectid   st_lengthshape   st_areashape      
 E09000001: 1   Barking and Dagenham: 1                :33   Min.   :294   Min.   : 8929   Min.   :  2897649  
 E09000002: 1   Barnet              : 1   Abertawe     : 0   1st Qu.:302   1st Qu.:28384   1st Qu.: 26797942  
 E09000003: 1   Bexley              : 1   Blaenau Gwent: 0   Median :310   Median :37664   Median : 37628571  
 E09000004: 1   Brent               : 1   Bro Morgannwg: 0   Mean   :310   Mean   :39255   Mean   : 47682317  
 E09000005: 1   Bromley             : 1   Caerdydd     : 0   3rd Qu.:318   3rd Qu.:46679   3rd Qu.: 56413925  
 E09000006: 1   Camden              : 1   Caerffili    : 0   Max.   :326   Max.   :74641   Max.   :150125298  
 (Other)  :27   (Other)             :27   (Other)      : 0                                                    

summary(BluePlaques)

Object of class SpatialPointsDataFrame
Coordinates:
             min      max
coords.x1 -0.477  0.21903
coords.x2  0.000 51.67830
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Number of points: 2812
Data attributes:
       id         
 Min.   :    1.0  
 1st Qu.:  711.8  
 Median : 6089.0  
 Mean   :10622.0  
 3rd Qu.:10358.2  
 Max.   :49190.0  
                  
                                                                                                                                                                              inscription  
 Frank Matcham (1854-1920) theatre architect designed this theatre                                                                                                                  :   3  
 Charlie Chester MBE                                                                                                                                                                :   2  
 Grim’s Dyke                                                                                                                                                                        :   2  
 Lillie Langtry 1852-1929 actress lived here                                                                                                                                        :   2  
  J. H. Greathead   Chief Engineer City and South London Railway. Inventor of the Travelling Shield that made possible the cutting of the tunnels of London's deep level tube system:   1  
  Joseph Wilson Swan and Sir James Dewar                                                                                                                                            :   1  
 (Other)                                                                                                                                                                            :2801  

Transform CRS from 4326 to 27700. (for Kernel Density Estimation, spatial data needs to be projected)

#now set up an EPSG string to help set the projection 
BNG = "+init=epsg:27700"
BluePlaques <- spTransform(BluePlaques, BNG)
BoroughMap <- spTransform(BoroughMap, BNG)

Plot Blue Plaques interactive map

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(BoroughMap) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(BluePlaques) +
  tm_dots(col = "blue")

Clip plaques to ‘LondonBorough’ boundaries

##first we'll remove any Plaques with the same grid reference as this will cause problems later on in the analysis...
BluePlaques <- remove.duplicates(BluePlaques)
#now just select the points inside London - thanks to Robin Lovelace for posting how to do this one, very useful!
BluePlaquesSub <- BluePlaques[BoroughMap,]
#check to see that they've been removed
tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(BoroughMap) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(BluePlaquesSub) +
  tm_dots(col = "blue")

Select one Borough for analysis

#OK, so just select the borough you are interested in. I'm going to try Haringey, you can choose any borough you like...
Borough <- BoroughMap[BoroughMap@data$lad15nm=="Haringey",]
#Check to see that the correct borough has been pulled out
tm_shape(Borough) +
  tm_polygons(col = NA, alpha = 0.5)

Subset Blue Plaques falling in Haringey

#clip the data to our single borough
BluePlaquesSub <- BluePlaques[Borough,]
#check that it's worked
tm_shape(Borough) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(BluePlaquesSub) +
  tm_dots(col = "blue")

We now have all of our data set up so that we can start the analysis using spatstat. The first thing we need to do is create an observation window for spatstat to carry out its analysis within - we’ll set this to the extent of the Haringey boundary

##now set a window as the borough boundary
window <- as.owin(Borough)
plot(window)

spatstat has its own set of spatial objects that it works with (one of the delights of R is that different packages are written by different people and many have developed their own data types) - it does not work directly with the SpatialPolygonsDataFrames, SpatialPointsDataFrames or sf objects that we are used to. For point pattern analysis, we need to create a point pattern (ppp) object.

#create a ppp object
BluePlaquesSub.ppp <- ppp(x=BluePlaquesSub@coords[,1],y=BluePlaquesSub@coords[,2],window=window)

Plot ppp object

plot(BluePlaquesSub.ppp,pch=16,cex=0.5, main="Blue Plaques Haringey")

Summarize point data - Kernel Density Estimation One way to summarise your point data is to plot the density of your points under a window called a ‘Kernel’. The size and shape of the Kernel affects the density pattern produced (more of this next week), but it is very easy to produce a KDE map from a ppp object using the density function.

plot(density(BluePlaquesSub.ppp, sigma = 500))

The sigma value sets the diameter of the Kernel (in the units your map is in - in this case, as we are in British National Grid the units are in metres). Try experimenting with different values of sigma to see how that affects the density estimate.

plot(density(BluePlaquesSub.ppp, sigma = 150))

1. Quadrat Analysis

The distribution of points in our study area differs from ‘complete spatial randomness’ - CSR? The most basic test of CSR is a quadrat analysis. We can carry out a simple quadrat analysis on our data using the quadrat count function in spatstat.

#First plot the points
plot(BluePlaquesSub.ppp,pch=16,cex=0.5, main="Blue Plaques in Haringey")
#now count the points in that fall in a 5 x 5 grid overlaid across the window
plot(quadratcount(BluePlaquesSub.ppp, nx = 5, ny = 5),add=T,col="red")

In our case here, want to know whether or not there is any kind of spatial patterning associated with the Blue Plaques in areas of London. If you recall from the lecture, this means comparing our observed distribution of points with a statistically likely (Complete Spatial Random) distibution, based on the Poisson distribution.

Using the same quadratcount function again (for the same sized grid) we can save the results into a table:

#run the quadrat count
Qcount<-data.frame(quadratcount(BluePlaquesSub.ppp, nx = 5, ny = 5))
#put the results into a data frame
QCountTable <- data.frame(table(Qcount$Freq, exclude=NULL))
#view the data frame
QCountTable

We don’t need the last row, so remove it

QCountTable <- QCountTable[-nrow(QCountTable),]
#check the data type in the first column - if it is factor, we will need to convert it to numeric
class(QCountTable[,1])

[1] "factor"

#oops, looks like it's a factor, so we need to convert it to numeric
vect<- as.numeric(levels(QCountTable[,1]))
vect <- vect[1:6]
QCountTable[,1] <- vect

OK, so we now have a frequency table - next we need to calculate our expected values. The formula for calculating expected probabilities based on the Poisson distribution is:

Pr(X=k)=λke/λk!

#calculate the total blue plaques (Var * Freq)
QCountTable$total <- QCountTable[,1]*QCountTable[,2]
#calculate mean
sums <- colSums(QCountTable[,-1])
sums

 Freq total 
   22    34 

Calculate Mean Poisson parameter

#and now calculate our mean Poisson parameter (lambda)
lambda <- sums[2]/sums[1]
#calculate expected using the Poisson formula from above - k is the number of blue plaques counted in a square and is found in the first column of our table...
QCountTable$Pr <- ((lambda^QCountTable[,1])*exp(-lambda))/factorial(QCountTable[,1])
#now calculate the expected counts and save them to the table
QCountTable$Expected <- round(QCountTable$Pr * sums[1],0)
QCountTable

Compare frequncy distributions

#Compare the frequency distributions of the observed and expected point patterns
plot(c(1,5),c(0,14), type="n", xlab="Number of Blue Plaques (Red=Observed, Blue=Expected)", ylab="Frequency of Occurances")
points(QCountTable$Freq, col="Red", type="o", lwd=3)

points(QCountTable$Expected, col="Blue", type="o", lwd=3)

Comparing the observed and expected frequencies for our quadrat counts, we can observe that they both have higher frequency counts at the lower end - something reminiscent of a Poisson distribution. This could indicate that for this particular set of quadrats, our pattern is close to Complete Spatial Randomness (i.e. no clustering or dispersal of points). But how do we confirm this?

To check for sure, we can use the quadrat.test function, built into spatstat. This uses a Chi Squared test to compare the observed and expected frequencies for each quadrat (rather than for quadrat bins, as we have just computed above). If the p-value of our Chi-Squared test is > 0.05, then we can reject a null hyphothesis that says “there is no complete spatial randomness in our data” (we will learn more about what a null-hypothesis is in a couple of weeks, but for the time being, just think about it as the opposite of a hypothesis that says our data exhibit complete spatial randomness). What we need to look for is a value for p > 0.05. If our p-value is > 0.05 then this indicates that we have CSR and there is no pattern in our points. If it is < 0.05, this indicates that we do have clustering in our points.

teststats <- quadrat.test(BluePlaquesSub.ppp, nx = 5, ny = 5)

Some expected counts are small; chi^2 approximation may be inaccurate

plot(BluePlaquesSub.ppp,pch=16,cex=0.1, main="Blue Plaques in Haringey")
plot(teststats, add=T, col = "red")

So we can see that the indications are there is no spatial patterning for Blue Plaques in Harrow - at least for this particular grid. Note the warning message - some of the observed counts are very small (0) and this may affect the accuracy of the quadrat test. Recall that the Poisson distribution only describes observed occurrances that are counted in integers - where our occurances = 0 (i.e. not observed), this can be an issue. We also know that there are various other problems that might affect our quadrat analysis, such as the modifiable areal unit problem.

In the new plot, we can see three figures for each quadrat. The top-left figure is the observed count of points; the top-right is the Poisson expected number of points; the bottom value is the Pearson residual value, or (Observed - Expected) / Sqrt(Expected).

2. Ripley’s K

One way of getting around the limitations of quadrat analysis is to compare the observed distribution of points with the Poisson random model for a whole range of different distance radii. This is what Ripley’s K function computes.
We can conduct a Ripley’s K test on our data very simply with the spatstat package using the kest function.

K <- Kest(BluePlaquesSub.ppp, correction="border")
plot(K)

The plot for K has a number of elements that are worth explaining. First, the Kpois(r) line in Red is the theoretical value of K for each distance window (r) under a Poisson assumption of Complete Spatial Randomness. The Black line is the estimated values of K accounting for the effects of the edge of the study area.

Where the value of K falls above the line, the data appear to be clustered at that distance. Where the value of K is below the line, the data are dispersed. “From the graph, we can see that up until distances of around 1300 metres, Blue Plaques appear to be clustered in Haringey, however, at around 1500 m, the distribution appears random and then dispersed between about 1600 and 2100 metres”.

Alternatives to Ripley’s K
There are a number of alternative measures of spatial clustering which can be computed in spatstat such as the G and the L functions - I won’t go into them now, but if you are interested, you should delve into the following references:

Bivand, R. S., Pebesma, E. J., & Gómez-Rubio, V. (2008). “Applied spatial data analysis with R.” New York: Springer.
Brundson, C., Comber, L., (2015) “An Introduction to R for Spatial Analysis & Mapping”. Sage.

https://research.csiro.au/software/wp-content/uploads/sites/6/2015/02/Rspatialcourse_CMIS_PDF-Standard.pdf

2. Density-based spatial clustering of applications with noise: DBSCAN

Quadrat and Ripley’s K analysis are useful exploratory techniques for telling us if we have spatial clusters present in our point data, but they are not able to tell us WHERE in our area of interest the clusters are occurring. To discover this we need to use alternative techniques. One popular technique for discovering clusters in space (be this physical space or variable space) is DBSCAN. For the complete overview of the DBSCAN algorithm, read the original paper by Ester et al. (1996) - http://www.aaai.org/Papers/KDD/1996/KDD96-037.pdf or consult the wikipedia page - https://en.wikipedia.org/wiki/DBSCAN

libraries

library(raster)

Attaching package: ‘raster’

The following objects are masked from ‘package:MASS’:

    area, select

The following objects are masked from ‘package:spatstat’:

    area, rotate, shift

The following object is masked from ‘package:nlme’:

    getData

library(fpc)

Error in library(fpc) : there is no package called ‘fpc’

DBSCAN requires you to input two parameters: 1. Epsilon - this is the radius within which the algorithm with search for clusters 2. MinPts - this is the minimum number of points that should be considered a cluster

“Based on the results of the Ripley’s K analysis earlier, we can see that we are getting clustering up to a radius of around 1200m, with the largest bulge in the graph at around 700m. Therefore, 700m is probably a good place to start and we will begin by searching for clusters of at least 4 points…”

So the DBSCAN analysis shows that for these values of eps and MinPts there are three clusters in the area I am analysing. Try varying eps and MinPts to see what difference it makes to the output.

No of course the plot above is a little basic and doesn’t look very aesthetically pleasing. As this is R and R is brilliant, we can always produce a much nicer plot by extracting the useful information from the DBSCAN output and use ggplot2 to produce a much cooler map…

rr library(ggplot2) #our new db object contains lots of info including the cluster each set of point coordinates belongs to, whether the point is a seed point or a border point etc. We can get a summary by just calling the object db

dbscan Pts=41 MinPts=3 eps=300
        0 1 2
border 29 1 1
seed    0 5 5
total  29 6 6

rr #if you open up the object in the environment window in RStudio, you will also see the various slots in the object, including cluster db$cluster

 [1] 0 0 1 1 0 1 0 2 2 0 2 1 0 2 0 0 0 0 2 0 0 1 0 0 0 0
[27] 0 0 2 0 0 0 0 0 0 0 0 1 0 0 0

rr #we can now add this cluster membership info back into our dataframe BluePlaquesSubPoints\(cluster <- db\)cluster #next we are going to create some convex hull polygons to wrap around the points in our clusters #use the ddply function in the plyr package to get the convex hull coordinates from the cluster groups in our dataframe chulls <- ddply(BluePlaquesSubPoints, .(cluster), function(df) df[chull(df\(coords.x1, df\)coords.x2), ]) # as 0 isn’t actually a cluster (it’s all points that aren’t in a cluster) drop it from the dataframe chulls <- subset(chulls, cluster>=1) #now create a ggplot2 object from our data dbplot <- ggplot(data=BluePlaquesSubPoints, aes(coords.x1,coords.x2, colour=cluster, fill=cluster)) #add the points in dbplot <- dbplot + geom_point() #now the convex hulls dbplot <- dbplot + geom_polygon(data = chulls, aes(coords.x1,coords.x2, group=cluster), alpha = 0.5) #now plot, setting the coordinates to scale correctly and as a black and white plot (just for the hell of it)… dbplot + theme_bw() + coord_equal()

2. Analysing Spatial Autocorrelation with Moran’s I, LISA and friends